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Abstract. Coherent techniques for searches of gravitational-wave bursts 
effectively combine data from several detectors, taking into account differences 



^•fj in their responses. The efforts are now focused on the maximum likelihood 



principle as the most natural way to combine data, which can also be used 

without prior knowledge of the signal. Recent studies however have shown that 

straightforward application of the maximum likelihood method to gravitational 

waves with unknown waveforms can lead to inconsistencies and unphysical results 

such as discontinuity in the residual functional, or divergence of the variance of 

the estimated waveforms for some locations in the sky. So far the solutions to 

_ . these problems have been based on rather different physical arguments. Following 

>— i^ ■ these investigations, we now find that all these inconsistencies stem from rank 

(_^ ' deficiency of the underlying network response matrix. In this paper we show 

nJ ' that the detection of gravitational-wave bursts with a network of interferometers 

C ^ ' belongs to the category of ill-posed problems. We then apply the method of 

^^3 1 Tikhonov regularization to resolve the rank deficiency and introduce a minimal 

C " 3 ■ regulator which yields a well-conditioned solution to the inverse problem for all 

""^"j ' locations on the sky. 

O . 

u 

{3J), PACS numbers: 04.80.Nn, 07.05.Kf, 95.55.Ym 



1. Introduction 

Efforts in searches for bursts with gravitational-wave detectors are now devoted 
to the development and testing of coherent data analysis techniques which do not 
require prior knowledge of the signal. Several such methods have been proposed in 
the past, beginning with the method of Giirsel and Tinto ^ which substantially 
predates all other techniques. The approach of Giirsel and Tinto is based upon 
explicit construction of a null data stream for two-polarization gravitational waves 
with completely arbitrary waveforms. Renewed interest in this method is motivated 
by the search for efficient vetoes of bursts of non-astrophysical origin |2]. Other 
methods which do not rely on the waveforms include coherent power filters [2j and 
cross-correlations of data from different detectors ^EJEI- Recent studies are however 
converging on the maximum likelihood method as the most natural way to combine 
data from a network of gravitational wave detectors. It was shown by Flanagan and 
Hughes that the maximum likelihood inference can, in principle, be used without 
any knowledge of the anticipated signal [7| . This observation led to the development 

I Current address: LIGO Hanford Observatory, PO Box 159, Richland, WA 99352, USA 



Tikhonov regularization for gravitational wave bursts 2 

of a data analysis technique known as excess power jS,. Subsequent extensions of 
this approach included Karhunen-Loeve expansions [HI and non-parametric adaptive 
filters ^U] ■ This groundwork helped building our confidence in the maximum likelihood 
method as a general framework for searches of gravitational- wave bursts with unknown 
waveforms. However, it was recently discovered that straightforward application of the 
maximum likelihood principle can lead to inconsistencies and unphysical results such 
as discontinuity in the residual functional called the two-detector paradox [111 IT^ , or 
unphysically large variations in the estimated signal-to- noise ratio ^H]- The proposed 
solutions included constraints and penalty functions derived from various physical 
arguments |12[I13| . Continuing these investigations, we now consider the problem from 
a very general point of view. In this paper, we show that all these inconsistencies and 
paradoxes arise because the inverse problem for bursts belongs to the category of ill- 
posed discrete (matrix) problems. In particular, its underlying matrix of coefficients 
suffers from rank deficiency. It is then natural to look for a solution within the 
Tikhonov regularization approach 1141 which is the general framework for solving ill- 
posed problems in mathematical physics. 

2. The inverse problem for bursts 

2.1. Network response 

Consider a network of m detectors located at different places on Earth. The response 
of the detectors to gravitational waves with two polarizations, /i+ and /ix , is given by 

Ut) ^ F,+{c^,e)h+{t) + F,^(c^,e)h^{i), (1) 

where Fi^ and Fi^ are the antenna-pattern functions, and (j) and 9 are the spherical 
angles of the source in the sky. In general, the data from the network contains both 
signal and noise (see figure^: 

X^{t)=i^{t + T,)+^,(t), (2) 

where the noise terms rji are assumed to be statistically independent among the 
detectors. The delays r^ = Ti{4),0) depend on the source location and are calculated 
with respect to a common reference, usually taken at the center of Earth. After 
changing the variables: t ~* t ~ Ti, we can write (|1I2|) in the matrix form: 

x{t\^,9)^A{4>,d)h{t)+r^{t\c^,d), (3) 

where x, h, and rj are column vectors and yl is ?n x 2 matrix, 

Fi+ Fix 
A= -. : , (4) 

which will be called the network response matrix. The notation x{t\(j), 6) implies that 
each component of the vector is shifted by its appropriate time delay: 

x,{t\(^,9)^x,[t^n{(t,,e)]. (5) 

For simplicity, we will often omit ((/>, 9) and write Q as 

x{t)=Ah{t)+r]{t). (6) 

Also, it will be convenient to view the matrix A as comprised of two vectors: 

^=[F+Fx], (7) 
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Figure 1. Simulated waveform injection. Left: a typical two-polarization 
waveform from numerical modeling of binary black hole coalescence ^]. Right: 
modeled detector outputs in Hl-L-G network for a source located at iji = —60° 
and 6» = 22°, with the matched-fiher SNR 4 of 8, 12 and 7. 



known as the range vectors of A J^ . Examples of numerical calculations in this paper 
correspond to LIGO Hanford (HI and H2), LIGO Livingston (L), VIRGO (V), and 
GEO-600 (G) detectors. 



2.2. Moore-Penrose inverse 

The inverse problem can be formulated as follows: given data from m detectors Xi{t) 
find the gravitational- wave amplitudes hi{t) and the source location in the sky (0,0) 
by solving 

A{<j>,e)h{t) = x{t\c^,e), (8) 

in which the data is contaminated with noise. In general, the problem cannot be solved 
exactly, and one looks for an approximate solution by minimizing the functional 

L[h]^\\x{t)-Ah{t)\\\ (9) 

where ||.|| stands for vector 2-norm: 



E 



■{t)dt, 



(10) 



defined over a suitable interval of observation T. The least squares (LSQ) functional 
L[h] is usually derived within the maximum likelihood approach. 

Variation of the functional L[h] with respect to hi{t) yields the normal equations, 
which in the matrix form can be written as 



Mh{t) ^ A^x{t), 
Then the solution is given by 

h{t) ^ A'^ x{t), 
where the 2 x m matrix 

At = M-^A^ 



where 



M = A^A. 



(11) 
(12) 
(13) 
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Figure 2. Typical search results for the simulated waveform injection shown 
in figure Left: injected (blue) and estimated (red) waveforms (after low-pass 
filtering). Right: sky map of the residual functional I17i . The position of the 
source (—60°, 22°) is marked by +. 



is known as the Moore-Penrose inverse or pseudoinverse of matrix A |17j . 
Consider now a situation when the data contains a signal, 

xit)=Aiq^s,Os)hsit)+rjit), (14) 

where {(j)s,9s) represent the source position in the sky. Then the solution acquires a 
non-zero expectation value: 

{h{t))^A^{cf>,0)A{^s,0s)hs{t). (15) 

At the true location of the source {h{t)) = hs{t), i.e. the solution given by the Moore- 
Penrose inverse yields an un-biased estimation of the gravitational- wave amplitudes. 
An example of a solution for h{t) obtained from noisy data is shown in figure |21 Note 
that H12|l provides the solution to the first part of the inverse problem: determination 
of h{t). In the next section we describe the solution to the second part of the problem: 
determination of the source position on the sky {(f>, 9). 

Introduction of the Moore-Penrose inverse allows us to define two complementary 
and orthogonal projector matrices: 

P = AA\ and Q = I - AA'' . (16) 

Note that PA — A and QA — which means that P projects onto the vector space of 
the range of matrix A and Q projects onto the null space of A"^ , which is a subspace 
complementary to the range space. 

2.3. The residual functional 

Substitution of the solution h{t) (|12|) into L[h] yields the residual functional: 

R^WQxf =x'^Qx, (17) 

where we used the fact that Q^ = Q. This functional will be used to find the location 
of the source in the sky. We will now prove that the expectation value of R{(j>, 0) 
reaches minimum at the true location of the source. Assume that the data contains a 
signal H14I) . then the residual functional is given by 

R^\\QA,hsf + 2ri^QA,h,+TfQi^, (18) 
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where Ag = A(0s,0s). The first term in H18|) contains only deterministic quantities 
and therefore its expectation value is the same as its current value. The expectation 
value of the second term in vanishes because the signal does not correlate with the 
noise. Finally, the expectation value of the third term is given by 

m „T 

{ri^Qri)= Y.'^-n Mt ~ n)i^,{t ~ T,)) dt 

= E Q^i I ^^Avf{t)) dt = a' tr(Q), (19) 

where we assumed, for simplicity, that the noise in the detectors is Gaussian and its 
variance is a. By transforming Q to the diagonal form, one can show that tr((5) = m—2. 
Therefore, 

{R{<j),e)) =\\QAsh,\\^+a^im-2). (20) 

At the true location of the source QAg ~ 0, and the function {R{(j>, 9)) reaches its 
minimum. This concludes the solution to the second part of the inverse problem: 
determination of the source position (see figure |2Jl. It is important to note that the 
function {R{4), 9)) may have several minima on the sky of which only one corresponds 
to the true location of the source. 

3. Examples of the inverse problem 

The simplest types of the inverse problem occur in networks of 2 and 3 detectors. 

3.1. m = 2 

In 2 dimensions the vectors F+ and F x span the entire vector space of columns of A, 
and there are no null vectors. In this case the inverse problem allows an exact solution 

h = A-^x, (21) 

for which the residual vanishes identically. It is worthwhile to obtain this result in a 
somewhat different way. Note that for m — 2 the matrix A is square and therefore 
At = A-i. Then 

P ^ AA'' = I, and Q = 0. (22) 

Consequently, R{4>, 9) = for every point on the sky, and no source localization is 
possible. In this case, every point in the sky, including the one which corresponds to 
the true source, can be viewed as a minimum of (i?((/), 9)). 

3.2. m == 3 

In 3 dimensions the vectors F+ and Fx span a 2-dimensional subspace of the vector 
space of columns of A. The complementary subspace, which is the space of null vectors, 
is therefore 1 dimensional, and is defined by one null vector, K. The null condition 
A^ K = can be written in vector notation as 

K-F+=K-Fx=0, (23) 

which implies that up to a multiplicative constant 

K = F+xFx. (24) 
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Knowing that Q projects onto the nuU space, we obtain 

n..-^^ and p.--5..-^SlISl (25) 

^^ ~ |Kp ' -^ ~ -^ |K|2 ■ ^ ' 

Then from the definition (|17() we find the residual functional, 

R^x^Qx=^^, (26) 

which is the minimization functional of Giirsel and Tinto. 

Next, we introduce two vectors to partition the Moore-Penrose inverse: 

^t = [H+HxF, (27) 

so that the solution for h l|12|) can be written as 

h+ = H+ • X, and h^ = Hx • x. (28) 

It is easy to show that the partition vectors are given by 

1 

-1 

m 

Indeed, from the definition of the Moore-Penrose inverse (|13|l we find that 

H+ = [A/-i]nr+ + [Af-i]i2Fx 

= (detM)-i [(Fx-Fx)F+-(F+-Fx)Fx] 

= (detM)-^ Fx X (F+ x Fx). (31) 

Combining this result with 

detM=|Kp, (32) 

we obtain H29|) . Similarly, one can derive (|30|l . Note that the solution for h, written 
in terms of H+ and Hx . is the waveform estimator of Giirsel and Tinto. 



H+ = PiFp Fx X K, (29) 

Hx = ^ F+ X K. (30) 



4. Difficulties witfi tlie LSQ minimization 

Direct minimization of the LSQ functional encounters various difficulties, some of 
which are briefly described in this section. 

4.I. Divergence of the expectation value of the solution 

Consider the solution given by the Moore-Penrose inverse (|12|l as a function of sky 
position: 

h{t) = Af-^(0, e)A^{(j>, 9) x{t). (33) 

As we have seen, on average, this solution correctly reproduces the waveform of the 
gravitational wave, provided that the estimated location of the source coincides with 
its true location. If, however, the estimated source location slightly deviates from the 
true location, the solution can be very different from the true waveform, especially 
near those places on the sky where the matrix M becomes singular. Further discussion 
of the variations of the estimated waveforms and solutions based on the penalty 
functional can be found in 13 . 
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4-2. Divergence of the variance of the solution 

Consider now the error in estimation of the gravitational- wave amplitudes which comes 
from the presence of noise in the data: 

5h{t)^A\4),e)r^{t). (34) 

Taking into account the fact that the noise in the detectors is uncorrelated, we obtain 

rri 

{5K{t)5h,{t')) = Y, [A%Mndvk{t)m{t')) 

^a^[M-%6{t-t'). (35) 

Therefore, the variance of the solution, 

T 

{ShUt))dt = a^[M-%, (36) 

/o 

diverges as the estimated source location approaches those places in the sky where the 

matrix M becomes singular. Solutions to this problem based on constraints applied 

to the waveforms can be found in [T^ . 

4-3. Divergence of the residual functional 

Finally, consider the residual functional 1)17(1 . According to this definition the residual 
diverges if matrix M becomes singular. This result becomes somewhat puzzling if 
we recall that in one particular case of singularity, namely when all detectors in the 
network are co-aligned, the residual in fact is well defined. Indeed, consider a network 
of m co-aligned detectors. In this case it is not possible to solve for both h+ and hx , 
and one can only solve for their linear combination: ^ ~ F+h+ + F^h^- Then the 
LSQ functional takes the form |H]: 

771 

i[e]-Eii^'W-^wii'' (37) 

with an obvious solution: £_ = — X^I^Li^*- I^ this case, the projectors P and Q can 
be defined as 

Pij = — , and Qij = Sij . (38) 

Consequently, the residual functional is given by 

R = x'^Qx= —\^ / [x^(t)- Xi (t)x, (t)] dt. (39) 

Mathematically, this is an altogether different inverse problem, and there is no obvious 
reason why the solution of 1(37(1 must be related to the solution of the original inverse 
problem 0. However, our physical intuition tells us that the network of nearly 
aligned detectors must behave very similarly to a network in which the same detectors 
become fully aligned. Yet, the residual calculated for nearly aligned detectors ((TTjl 
does not approach the one calculated for co-aligned detectors (139(1 in the limit when 
the detectors become perfectly aligned. This discontinuity of the residual functional 
becomes even more striking in two-detector networks for which R = (see section ITT|l 
for all detector orientations, no matter how close to perfect alignment they are. In 
this case the discrepancy was called the two- detector paradox ;11..12.il3| . 
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Figure 3. Condition number as a function of sky position (colormap in log scale). 
Left: LIGO-only network (H1-H2-L), Right: LIGO-GEO network (H1-H2-L-G). 
Note that close alignment of LIGO interferometers leads to very high condition 
number for some locations on the sky (> 10^), and the inclusion of GEO detector 
results in significant reduction of the condition number (> 10^). 



5. Rank deficiency 

All the above problems originate from rank deficiency of the network response matrix 
A. Recall that the nominal rank of A is 2. The rank of A drops to 1 if all the rows of 
A become proportional to each other, which is known as full row degeneracy. Column 
degeneracy occurs when the two columns of matrix A become proportional to each 
other. Note that full row degeneracy implies column degeneracy of A and vice versa. 
In other words, either degeneracy results in coUinearity of the range vectors: 



/5Fh 



(40) 



where (i is real. If this condition is satisfied, the rank of A is 1. The simplest case of 
rank deficiency occurs when all detectors in the network are co-aligned and therefore 
all the rows of matrix A are equal. 

In practice, we are seldom concerned with the exact coUinearity of F+ and Fx as 
the problem already occurs when these vectors are close to being coUinear. This can 
be seen from the fact that 



det M 



vriFxr-(F+-Fx)"^o, (41) 

as the two range vectors approach coUinearity, making the Moore-Penrose inverse 
divergent. Quantitatively, the degree to which the inversion of A becomes ill defined 
is described by the condition number |18| : 

cond{A) = \\A\\-\\A^l (42) 

where ||.|| stands for matrix 2-norm.§ Perfectly-invertible orthogonal matrices have 
condition number of 1. Large condition numbers indicate ill-defined inversion. Figure|31 
shows condition numbers for two examples of detector networks. The locations in the 
sky with large condition number correspond to rank-deficient network response matrix. 

§ The 2-norm of matrix A is defined as the maximum value of the 2-norm of vector Av under 
condition that |jv|| = 1. 
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6. Tikhonov regularization 

Several techniques are available in applied mathematics to address rank deficiency 
of the coefScient matrix in the LSQ problem ^Sl- One of the most commonly used 
and best understood techniques is the Tikhonov regularization method ^3]. The key 
idea in this method is to introduce a regularization functional (regulator) f2[/i(i)] with 
parameter 5 > such that the modified LSQ functional, 

Lg[h] = \\x{t) -Ah{t)f + g n[h], (43) 

no longer suffers from rank deficiency. Consider for example a quadratic regulator, 

n[h] = h'^nh= Y^ / n^Jh^{t)hJ{t)dt, (44) 

where flij is a symmetric 2x2 matrix. Quadratic regulators preserve the linearity 
of the inverse problem and therefore have an advantage over other forms. Then the 
solution of the inverse problem becomes 

h = M-^A^x, where Mg = M + gfl. (45) 

Therefore, the regularized version of the Moore-Penrose inverse is 

Al = M-'A'^, (46) 

which leads to the following generalization of matrices P and Q: 

Pg=AAl, and Qg = I - AAl (47) 

In general, these matrices no longer satisfy the property of projectors. However, the 
connection between the matrix Q and the residual functional still holds. Indeed, 
substituting the solution (|I5|) into the modified LSQ functional (|^ . we obtain 

Rg=x^QgX, (48) 

which is equivalent to (|17|l despite the presence of a regulator in Lg [h] . 

The role of SI is to render the inverse of Mg well defined when M is nearly singular. 
This resolves the problems associated with singularities of A/, described in sections 
14.11 and 14.21 Introduction of the regulator also solves the problem of discontinuity of 
the residual, described in section lTTSl Indeed, for .g 7^ the residual Rg is a continuous 
function of the detector alignment, and no divergence occurs in Rg when the matrix 
A becomes rank deficient (|40|) . Furthermore, explicit calculations show that in this 
case 

^^'^'' " a + gdetn ^'^'' ^^^^ 

where f is a unit vector along F+, and 

a=(/32r!n-2/3f)i2 + r!22)|F+p. (50) 

If the detectors in the network become co-aligned, f,- = -k= and therefore 

a 1 

^^'^'' = a + gdetn m' ^^^' 

One can see from this expression that Pg reduces to P in (|38|l . in the limit g ^ 0, and 
consequently, Qg reduces to Q in H38|) . Hence, the residual functional H48|l reduces to 
that in ^. 
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Figure 4. Condition number for LIGO-VIRGO networ]^ {H1-H2-L-V) as a 
function of sl<y position for different values of g (colormap in log scale). Complete 
regularization takes place for 3=1 (bottom-right) in which case the condition 
number is 1. 



7. Optimum condition regulator 

As we have seen, the degree to which matrix A becomes non-invertible, measured by 
the condition number, strongly depends on the sky location. It is therefore desirable 
to construct a regulator which is a function of sky position. Particularly useful will be 
regulators which adjust themselves to higher condition number, always guaranteeing 
well-defined inversion of A. Here we construct one example of such a regulator. 
Consider matrix M in the space of its eigenvectors: 





M,, 







M2 



(52) 



It is easy to show that its eigenvalues, /ii and /Lt2, are always positive and one is always 
greater than the other, e.g. /ii > /i2. Since the purpose of a regulator is to prevent 
singularities which occur when /ii^2 ~^ 0, it is sufficient to consider U, which is diagonal 
in the space of eigenvectors of A'/: 

wi 

072 

where we assume that cJi.2 > so that U,[h] is positive definite. 



n. 



(53) 
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Figure 5. The residual functional EH}, for LIGO-VIRGO network (H1-H2- 
L-V) with no regularization (left) and full regularization (right). The injected 
waveforms (figure 1) are amplified to produce the signals in detectors with the 
SNR of 8, 8, 12, and 10 for source location (-60°, 22°) (marked by +). Note a 
significant reduction in the degeneracy (blue area) introduced by the regulator. 



In most cases of interest only one of the eigenvalues, /i2, can be singular. If both 
eigenvalues vanish then tr(Af ) = 0. This would imply that |F-|.p -t- |Fx P = 0, which, 
in turn, implies that all components of vectors F+ and Fx vanish. In a given detector 
both antenna-pattern functions vanish only if the source is located in the plane of 
the detector, on the bisector of the detector arms, or on the normal to the bisector. 
Apart from the two Hanford interferometers, a general network of detectors does not 
have common arm bisectors or bisector normals, which is why both eigenvalues cannot 
vanish simultaneously. 

A singularity in which one of the eigenvalues approaches zero whereas the other 
remains finite is known as a finite gap |18j . In this case, it is sufhcient to regularize 
only the smallest eigenvalue /i2- We therefore limit ourselves to regulators in which 
uji = 0, and which are often called semi-norms in the space of solutions. Note that the 
parameter g becomes redundant; it can be absorbed into lj2- Nevertheless, we retain 
g so that we can control the strength of the regulator, and assume that the maximum 
value of g is 1. We can now construct a regulator which makes the matrix A fully 
invertible over the entire sky. Quantitatively, this means that for .g = 1 

1411 = 1, 



for all 6 and 



cond{A) = \\A\\ ,,.g, 
Taking into account that 



l^ll 



Ml, 



and 



|At|i=max 



1 



/il ^2 + W2 



we find that the minimum value for such LO2 is given by 

M2- 



(54) 



(55) 



(56) 



t^2 ^ \/MiM2 

In other words, any regulator of the semi-norm type in which 0^2 > 0J2 yields the 
inversion of A with condition number of 1. 

Figure 01 shows the improvement in the condition number for the LIGO-VIRGO 
network which results from this regularization. This will imply improvement in the 
accuracy of the solution for hit). Note that the introduction of the regularization also 
improves localization of the source on the sky, as shown in figure [5| Further analysis of 
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the role of regularization will be given elsewhere. We conclude with the reminder that 
regularization, by its nature, introduces a bias and therefore the optimal approach 
must be a trade-off between the bias and the error due to noise. 

8. Conclusion 

We have shown that detection of gravitational-wave bursts of unknown waveforms is 
a linear inverse problem which becomes ill posed because of the rank deficiency of 
the underlying network response matrix. Following the general scheme of Tikhonov 
regularization, we introduced a semi-norm regulator in the space of solutions which 
guarantees full inversion of the network response matrix over the entire sky. The 
analysis presented here is general and applies to any network of interferometric 
gravitational- wave detectors. The problem of rank deficiency is particularly important 
from a practical point of view because the condition number for the LIGO-only network 
can be extremely high (> 10^) and regularization of the response matrix will play a 
crucial role in stabilizing the solution of the inverse problem. For networks in which 
LIGO detectors are accompanied by VIRGO or GEO interferometer the condition 
number is substantially better (> 10^) and yet still in need of regularization. 
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